When to use a time series vs. a linear regression?
If you have continuous target variable, then it is a regression problem. For instance, in flight trials we have the flight distance to predict, which is continuous. Hence this becomes a regression problem. About time series, when the datapoints are time dependent, then it becomes a time series problem. Each data point has an order and is, typically, related to the data points before and after by some underlying process.
In turn, a times series regression is a statistical method for predicting a future response based on the response history (known as autoregressive dynamics) and the transfer of dynamics from relevant predictors. Time series regression is commonly used for modeling and forecasting of economic, financial, and biological systems.
There are three concepts then to keep in mind:
Stationary: if any trends or patterns in the time series are independent with time. (A linear regression equivalent to a time series at stationarity). There is strong and weak stationary.
Non-Stationary: trends dependent in time in the mean or statistical deviation of the data points.
Autocorrelation: “memory”; the degree to which time series values in period(t) are related to time series values in periods (t+1, t+2, t+3…). For this, we can test to see how long certain values last in a system.
There are also three words and their definitions to remember:
Questions about time series often asked include,
So how can you answer these questions through a time series regression and code? And where do you begin to test stationarity or visualize autocorrelation? And finally, how do you handle non-stationarity?
source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat_working_on/Rsrc/"
script_names = c("compare_models.R","regression_output.R", "clean_morph_data.R",
"AICprobabilities.R")
for (script in script_names) {
path = paste0(source_path, script)
source(path)
}
Extract date information in order to prep it for converting it into a datetime object read by the xts() function.
Merge close dates together in order to have as equally spaced datapoints as possible. The acf() assumes that the data are regular spaced. The acf computes estimates of the autocovariance or autocorrelation function. We need this because it will tell us, if there are patterns, what the mathematical function of those patterns are.
data_list <- read_morph_data("data/allmorphology9.21.20.csv")
raw_data = data_list[[1]]
all_bugs = nrow(raw_data)
# data_long = data_list[[2]] # need to fix this
# Remove individuals with torn wings first
raw_data$drop <- FALSE
for(row in 1:nrow(raw_data)){
if(length(unlist(strsplit(strsplit(paste("test ", raw_data$notes[row], " test", sep=""), "torn")[[1]], "wing")))>2){
#browser()
raw_data$drop[row] <- TRUE
}
}
raw_data <- raw_data[raw_data$drop==FALSE,]
clean_bugs = nrow(raw_data)
cat("number of bugs with torn wings:", all_bugs - clean_bugs, "\n\n")
## number of bugs with torn wings: 141
# Datetime
raw_data$date <- paste(raw_data$month, raw_data$year, sep="/")
raw_data$datetime <- as.yearmon(raw_data$date, "%B/%Y")
raw_data$datetime <- as.factor(raw_data$datetime)
n_missing_dates = nrow(raw_data[is.na(raw_data$datetime),])
# merge May 2015 with April 2015 because very few bugs were collected in May 2015.
# then merge April 2013 with May 2013 to make the time datapoints more evenly distanced
raw_data$date[raw_data$date == "May/2015"] = "April/2015"
raw_data$date[raw_data$date == "May/2013"] = "April/2013"
# convert to yearmon object and then factor
raw_data$datetime <- as.yearmon(raw_data$date, "%B/%Y")
raw_data$datetime <- as.factor(raw_data$datetime)
cat("number of missing dates:", n_missing_dates, "\n\n")
## number of missing dates: 296
unique(raw_data$datetime)
## [1] Apr 2013 Apr 2014 Apr 2015 Dec 2013 Aug 2017 Dec 2016 Sep 2018 May 2019
## [9] Oct 2019 Feb 2020 <NA>
## 10 Levels: Apr 2013 Dec 2013 Apr 2014 Apr 2015 Dec 2016 Aug 2017 ... Feb 2020
Looks like we have a measurement at least once per year, which is good.
Remove missing wing and datetime values.
Compute the wing length averages and generate datetime objects using as.Date(). In order to be processed by the as.Date() function, all time objects need a date, month, and year.
Use the xts & zoo R Libraries to read cleanly index the data by a formal time object (collection_time).
Optional: include major events that occurred in the 2010’s. This could be events that you find biologically significant to your questions. E.g. how did major hurricanes in Florida impact average soapberry bug wing length across Florida?
Hurricane Source: https://en.wikipedia.org/wiki/List_of_Florida_hurricanes#2000–present
# remove NA dates
d = raw_data %>%
filter(!is.na(wing), !is.na(datetime))
# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
# events
FL_major_hurr = c("Sep 2017 10", "Oct 2018 10") # Irma (s), Michael (n)
hurr_dates = as.Date(FL_major_hurr, "%b %Y %d")
events <- xts(c("Irma (S. FL)", "Michael (N. Fl)"),
hurr_dates)
Create xts-zoo object and plot the time series data. Add any events to the plot with addEventLines().
Calculate the Augmented Dickey-Fuller Test to test for stationarity where the null hypothesis is non-stationarity and the alternative hypothesis is stationarity.
Finally, use an ACF (autocorrelation function) and PCF (partial-autocorrelation function) plot to identify temporal dependence in the data. Autocorrelation measures the linear relaionship between lagged values of a time series.
\(R_s = Corr(x_t, x_{t+s})\) for lag \(s\).
There are two horizonal, blue, dashed lines as well. Those represent the significance threshold, where only the spikes that exceed this dashed line are considered significant. In other words, these plots describe how well the present value of the series is related with its past values. Since, none do significantly, then there is no AR (a present value of the time series cannot be obtained using previous values of the same time series).
check_stationarity = function(dx) {
# dx is your xts zoo time series object
plot(dx)
#addEventLines(events, pos=2, srt=90, col="red")
par(mfrow=c(1,2))
print(adf.test(dx[,1]))
Acf(dx[,1], main='') # ACF
Acf(dx[,1], type="partial", main='') #PACF
par(mfrow=c(1,1))
}
##
## Augmented Dickey-Fuller Test
##
## data: dx[, 1]
## Dickey-Fuller = -2.5124, Lag order = 2, p-value = 0.3772
## alternative hypothesis: stationary
Interpretation:
Stationarity can be defined in precise mathematical terms, but for our purpose we mean a flat looking series, without trend, constant variance over time, a constant autocorrelation structure over time and no periodic fluctuations (seasonality). So the high p value means we do have non-stationarity. Similarly, the augmented Dickey–Fuller (ADF) statistic, used in the test, is a negative number. The more negative it is, the stronger the rejection of the null hypothesis. This is not negative enough.
Finally, the lag length (lag order) is how many terms back down the AR process you want to test for serial correlation. This is a non-stationary with AR(2).
Differencing is when a new series is constructed where the value at the current time step is calculated as the difference between the original observation and the observation at the previous time step.
detrend = function(dx) {
dx$diff <- diff(dx[,1]) # "regressing out"
dx <- dx[-1,]
plot(dx$diff)
par(mfrow=c(1,2))
print(adf.test(dx$diff))
Acf(dx$diff, main='') # ACF
Acf(dx$diff, type="partial", main='') #PACF
par(mfrow=c(1,1))
}
dedrift = function(dx) {
dx$logv <- log(dx[,1])
dx <- dx[-1,]
plot(dx$logv)
par(mfrow=c(1,2))
print(adf.test(dx$logv))
Acf(dx$logv, main='') # ACT
Acf(dx$logv, type="partial", main='') #PACF
par(mfrow=c(1,1))
}
dedriftrend = function(dx) {
dx$logdiff <- diff(log(dx[,1]))
dx <- dx[-1,]
plot(dx$logdiff)
par(mfrow=c(1,2))
print(adf.test(dx$logdiff))
Acf(dx$logdiff, main='') # ACT
Acf(dx$logdiff, type="partial", main='') #PACF
par(mfrow=c(1,1))
}
detrend(wing_mm) # this is stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$diff
## Dickey-Fuller = -3.7455, Lag order = 2, p-value = 0.03961
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
Interpretation: This time series is stationary! And as a reminder, a data that is just noise will be stationary.
dedrift(wing_mm) # this is stationary
## Warning in adf.test(dx$logv): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -6.4449, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedriftrend(wing_mm) # this is stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -4.2683, Lag order = 2, p-value = 0.0139
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
# Get only bugs with long wings
data_long<-raw_data[raw_data$w_morph=="L",]
# Calculate wing2body ratio for bugs with long wings
data_long$wing2body <- data_long$wing/as.numeric(data_long$body)
# remove NA dates
d = data_long %>%
filter(!is.na(wing2body), !is.na(datetime))
# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
ratio = xts(ratio_avg, date)
colnames(ratio) <- "wing2body"
function() {plot(ratio)}
## function() {plot(ratio)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(ratio$wing2body) # this is not stationary
## Warning in adf.test(ratio$wing2body): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ratio$wing2body
## Dickey-Fuller = 10.559, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(ratio$wing2body, main='') # ACT
Acf(ratio$wing2body, type="partial", main='') #PACF
detrend(ratio) # this is stationary
## Warning in adf.test(dx$diff): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$diff
## Dickey-Fuller = -7.7839, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedrift(ratio) # this is not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -1.0499, Lag order = 1, p-value = 0.9133
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedriftrend(ratio)
## Warning in adf.test(dx$logdiff): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -7.6591, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1
# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
filter(!is.na(wing_morph_binom), !is.na(datetime))
# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
morph = xts(freq_avg, date)
colnames(morph) <- "wing_morph_freq"
function() {plot(morph)}
## function() {plot(morph)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(morph$wing_morph_freq) # this is not stationarity
##
## Augmented Dickey-Fuller Test
##
## data: morph$wing_morph_freq
## Dickey-Fuller = -1.8526, Lag order = 2, p-value = 0.6285
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(morph$wing_morph_freq, main='') # ACT
Acf(morph$wing_morph_freq, type="partial", main='') #PACF
detrend(morph) # this is not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$diff
## Dickey-Fuller = -1.4413, Lag order = 2, p-value = 0.7852
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedrift(morph) # this is stationary
## Warning in adf.test(dx$logv): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -27.979, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedriftrend(morph) # this is not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -2.0047, Lag order = 2, p-value = 0.5706
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
females = raw_data[raw_data$sex=="F",]
males = raw_data[raw_data$sex=="M",]
# remove NA dates
d = females %>%
filter(!is.na(wing), !is.na(datetime))
# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
## function() {plot(wing_mmf)}
##
## Augmented Dickey-Fuller Test
##
## data: wing_mmf$wing
## Dickey-Fuller = -2.8815, Lag order = 2, p-value = 0.2366
## alternative hypothesis: stationary
Looks like there’s not a linear trend, but changes in variation over time.
dedrift(wing_mmf)
## Warning in adf.test(dx$logv): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -5.4041, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
# remove NA dates
d = data_long %>%
filter(!is.na(wing2body), !is.na(datetime))
d = d[d$sex=="F",]
# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
ratiof = xts(ratio_avg, date)
colnames(ratiof) <- "wing2body"
function() {plot(ratiof)}
## function() {plot(ratiof)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(ratiof$wing2body) # this is not stationarity
## Warning in adf.test(ratiof$wing2body): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ratiof$wing2body
## Dickey-Fuller = 1.4005, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(ratiof$wing2body, main='') # ACT
Acf(ratiof$wing2body, type="partial", main='') #PACF
dedrift(ratiof) # still not stationary
## Warning in adf.test(dx$logv): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = 0.051299, Lag order = 1, p-value = 0.99
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
dedriftrend(ratiof)
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -2.4151, Lag order = 1, p-value = 0.4143
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red") # *still* not stationary
raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1
# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
filter(!is.na(wing_morph_binom), !is.na(datetime))
d = d[d$sex=="F",]
# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
morphf = xts(freq_avg, date)
colnames(morphf) <- "wing_morph_freq"
function() {plot(morphf)}
## function() {plot(morphf)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(morphf$wing_morph_freq) # this is not stationary
##
## Augmented Dickey-Fuller Test
##
## data: morphf$wing_morph_freq
## Dickey-Fuller = -0.77584, Lag order = 2, p-value = 0.952
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(morphf$wing_morph_freq, main='') # ACT
Acf(morphf$wing_morph_freq, type="partial", main='') #PACF
detrend(morphf) # still not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$diff
## Dickey-Fuller = -0.69068, Lag order = 2, p-value = 0.9591
## alternative hypothesis: stationary
dedriftrend(morphf) # still not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -0.83913, Lag order = 2, p-value = 0.9442
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
# remove NA dates
d = males %>%
filter(!is.na(wing), !is.na(datetime))
# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
## function() {plot(wing_mmm)}
##
## Augmented Dickey-Fuller Test
##
## data: wing_mmm$wing
## Dickey-Fuller = -0.83612, Lag order = 2, p-value = 0.9447
## alternative hypothesis: stationary
dedrift(wing_mmm) # not stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -1.4491, Lag order = 2, p-value = 0.7822
## alternative hypothesis: stationary
dedriftrend(wing_mm) # now stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logdiff
## Dickey-Fuller = -4.2683, Lag order = 2, p-value = 0.0139
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
# remove NA dates
d = data_long %>%
filter(!is.na(wing2body), !is.na(datetime))
d = d[d$sex=="M",]
# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
ratiom = xts(ratio_avg, date)
colnames(ratiom) <- "wing2body"
function() {plot(ratiom)}
## function() {plot(ratiom)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(ratiom$wing2body) # this is not stationarity
##
## Augmented Dickey-Fuller Test
##
## data: ratiom$wing2body
## Dickey-Fuller = -0.44044, Lag order = 2, p-value = 0.9776
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(ratiom$wing2body, main='') # ACT
Acf(ratiom$wing2body, type="partial", main='') #PACF
detrend(ratiom) # this is stationary
## Warning in adf.test(dx$diff): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dx$diff
## Dickey-Fuller = -7.954, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1
# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
filter(!is.na(wing_morph_binom), !is.na(datetime))
d = d[d$sex=="M",]
# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
morphm = xts(freq_avg, date)
colnames(morphm) <- "wing_morph_freq"
function() {plot(morphm)}
## function() {plot(morphm)}
addEventLines(events, pos=2, srt=90, col="red")
adf.test(morphm$wing_morph_freq) # this is not stationary
##
## Augmented Dickey-Fuller Test
##
## data: morphm$wing_morph_freq
## Dickey-Fuller = -1.3425, Lag order = 2, p-value = 0.8229
## alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(morphm$wing_morph_freq, main='') # ACT
Acf(morphm$wing_morph_freq, type="partial", main='') #PACF
dedrift(morphm) # close! But marginally stationary
##
## Augmented Dickey-Fuller Test
##
## data: dx$logv
## Dickey-Fuller = -3.5464, Lag order = 2, p-value = 0.05744
## alternative hypothesis: stationary
addEventLines(events, pos=2, srt=90, col="red")
# remove NA dates
d = males %>%
filter(!is.na(wing), !is.na(datetime))
# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)
# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")
## function() {plot(wing_mmm)}
##
## Augmented Dickey-Fuller Test
##
## data: wing_mmm$wing
## Dickey-Fuller = -0.83612, Lag order = 2, p-value = 0.9447
## alternative hypothesis: stationary